The previous chapter described the integration of two pathways. We will now add yet another set of reactions to the previously integrated pathways. The AMP input and output that was described earlier by two simple reactions across the system boundary represent a network of reactions that are called purine degradation and salvage pathways. The salvage pathways recycle a pentose through the attachment of an imported purine base to effectively reverse AMP degradation. AMP degradation and salvage are low flux pathways but have several genetic mutations in the human population that cause serious diseases. Thus, even though fluxes through these pathways are low, their function is critical. As we have seen, many of the long-term dynamic responses are determined by the adjustment of the total adenylate cofactor pool. In this chapter, we integrate the degradation and salvage pathways to form a core metabolic network for the human red blood cell.
MASSpy will be used to demonstrate some of the topics in this chapter.
from mass import (
MassModel, MassMetabolite, MassReaction,
Simulation, MassSolution, strip_time)
from mass.io.json import load_json_model
from mass.util.matrix import nullspace, left_nullspace, matrix_rank
from mass.visualization import (
plot_time_profile, plot_phase_portrait, plot_tiled_phase_portraits,
comparison_plot)
Other useful packages are also imported at this time.
from os import path
from cobra import DictList
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import sympy as sym
Some options and variables used throughout the notebook are also declared here.
pd.set_option("display.max_rows", 100)
pd.set_option("display.max_columns", 100)
pd.set_option('display.max_colwidth', -1)
pd.options.display.float_format = '{:,.3f}'.format
S_FONT = {"size": "small"}
L_FONT = {"size": "large"}
INF = float("inf")
AMP metabolism forms a sub-network in the overall red blood cell metabolic network. We will first consider its properties before we integrate it with the combined glycolytic and pentose pathway model.
AMP is simultaneously degraded and synthesized, creating a dynamic balance (Figure 12.1). The metabolites that are found in these pathways and that are over and above those that are in the integrated glycolytic and pentose pathway network are listed in Table 12.1. To build this sub-network, prior to integration with the glycolytic and pentose pathway network, we need consider a sub-network comprised of the metabolites shown in Table 12.1. The reactions that take place in the degradation and biosynthesis of AMP are shown in Table 12.2.
Figure 12.1: The nucleotide metabolism considered in this chapter. There are two biochemical ways in which AMP is degraded and two ways in which it is synthesized.
filepath = path.realpath(path.join(
"models", "SB2_AMPSalvageNetwork.json"))
ampsn = load_json_model(filepath)
AMP can be degraded either by dephosphorylation or by deamination. In the former case, adenosine (ADO) is formed and can cross the cell membrane to enter plasma. In the latter case, IMP is formed and can then subsequently be dephosphorylated to form inosine (INO). INO can then cross the cell membrane, or be further degraded to form a pentose (R1P) and the purine base hypoxanthine (HYP), the latter of which can be exchanged with plasma. HYP will become uric acid, the accumulation of which causes hyperuricemia (gout).
The red blood cell does not have the capacity for de novo adenine synthesis. It can synthesize AMP in two different ways. First, it can phosphorylate ADO directly to form AMP using an ATP to ADP conversion. This ATP use creates an energy load met by glycolysis. A second process to offset degradation of adenine is the salvage pathway, in which adenine (ADE) is picked up from plasma and is combined with phosphoribosyl diphosphate (PRPP) to form AMP. PRPP is formed from R5P using two ATP equivalents, where the R5P is formed from isomerization of R1P that is formed during INO degradation. In this way, the pentose is recycled to re-synthesize AMP at the cost of two ATP molecules. The R5P can also come from the pentose phosphate pathway in an integrated network.
Table 12.1: Intermediates of AMP degradation and synthesis, their abbreviations and steady state concentrations. The concentrations given are those typical for the human red blood cell. The index on the compounds is added to that for the combined glycolysis and pentose phosphate pathway, (Table 11.1).
Table 12.2: AMP degradation and biosynthetic pathway enzymes and transporters, their abbreviations and chemical reactions.The index on the reactions is added to that for the combined glycolysis and pentose pathways (Table 11.2)
Table 12.3: The elemental composition and charges of the AMP Salvage Network intermediates. This table represents the matrix $\textbf{E}.$
Table 12.4: The elemental and charge balance test on the reactions. All internal reactions are balanced. Exchange reactions are not.
for boundary in ampsn.boundary:
print(boundary)
Five pathway vectors chosen to span the null space are shown in Table 12.5. They can be divided into groups of degradative and biosynthetic pathways.
Table 12.5: The calculated null space pathway vectors for the stoichiometric matrix for the AMP Salvage Network.
The first three pathways, $\textbf{p}_1$, $\textbf{p}_2$, and $\textbf{p}_3$ are degradation pathways of AMP. The first pathway degrades AMP to ADO via dephosphorylation. The second pathway degrades AMP to IMP via AMP deaminase followed by dephosphorylation and secretion of INO. The third pathway degrades AMP first with dephosphorylation to ADO followed by deamination to INO and its secretion. These are shown graphically on the reaction map in Figure 12.2. Note that the second and third pathways are equivalent overall, but take an alternative route through the network. These two pathways are said to have an equivalent input/output signature (Systems Biology: Properties of Reconstructed Networks).
The fourth pathway, $\textbf{p}_4$, shows the import of adenosine (ADO) and its phosphorylation via direct use of ATP. Thus, the cost of the AMP molecule now relative to the plasma environment is one high-energy phosphate bond. In the previous chapters, AMP was directly imported and did not cost any high-energy bonds in the defined system. Pathway four is essentially the opposite of pathway one, except it requires ATP for fuel. If one sums these two pathways, a futile cycle results.
The fifth pathway, $\textbf{p}_5$, is a salvage pathway. It takes the INO produced through degradation of AMP, cleaves the pentose off to form and secrete hypoxanthine (HYP), and recycles the pentose through the formation of PRPP at the cost of two ATP molecules. PPRP is then combined with an imported adenine base to form AMP. This pathway is energy requiring, needing two ATPs. As we detail below, some of the chemical reactions of this pathway change when it is integrated with glycolysis and the pentose pathway, as we have modified the PRPP synthase reaction to make it easier to analyze this sub-network.
Figure 12.2: The graphical depiction of the five chosen pathway vectors of the null space of the stoichiometric matrix, shown in Table 12.1.5. From left to right: The first three pathways $\textbf{p}_1$, $\textbf{p}_2$, and $\textbf{p}_3$, are degradation pathways, while the $\textbf{p}_4$ is a direct synthesis pathway and $\textbf{p}_5$ is a salvage pathway.
The left null space is one dimensional. It has one time invariant: ATP + ADP. In this sub-network, this sum acts like a conserved cofactor.
Table 12.6: The left null space composed of the time invariant ATP + ADP pool.
All of the properties of the stoichiometric matrix can be conveniently summarized in a tabular format, Table 12.7. The table succinctly summarizes the chemical and topological properties of $\textbf{S}$. The matrix has dimensions of 15x19 and a rank of 14. It thus has a 5 dimensional null space and a 1 dimensional left null space.
Table 12.7: The stoichiometric matrix for the AMP Salvage Network seen in Figure 12.1. The matrix is partitioned to show the intermediates (yellow) separate from the cofactors and to separate the exchange reactions and cofactor loads (orange). The connectivities, $\rho_i$ (red), for a compound, and the participation number, $\pi_j$ (cyan), for a reaction are shown. The second block in the table is the product $\textbf{ES}$ (blue) to evaluate elemental balancing status of the reactions. All exchange reactions have a participation number of unity and are thus not elementally balanced. The last block in the table has the five pathway vectors (purple) for the AMP Salvage Network. These vectors are graphically shown in Figure 12.2. Furthest to the right, we display the time invariant pools (green) that span the left null space.
The null space is five dimensional. We thus have to specify five fluxes to set the steady state.
In a steady state, the synthesis of AMP is balanced by degradation, that is $v_{SK_{amp}}=0$. Thus, the sum of the flux through the first three pathways must be balanced by the fourth to make the AMP exchange rate zero. Note that the fifth pathway has no net AMP exchange rate.
The fifth pathway is uniquely defined by either the exchange rate of hypoxanthine or adenine. These two exchange rates are not independent. The uptake rate of adenine is approximately 0.014 mM/hr (Joshi, 1990).
The exchange rate of adenosine would specify the relative rate of pathways one and four. The rate of $v_{ADNK1}$ is set to 0.12 mM/hr, specifying the flux through $\textbf{p}_{4}$. The net uptake rate of adenosine is set at 0.01 mM/hr, specifying the flux of $\textbf{p}_{1}$ to be 0.11 mM/hr.
Since $\textbf{p}_{1}$ and $\textbf{p}_{4}$ differ by 0.01 mM/hr in favor of AMP synthesis, it means that the sum of $\textbf{p}_{2}$ and $\textbf{p}_{3}$ has to be 0.01 mM/hr. To specify the contributions to that sum of the two pathways, we would have to know one of the internal rates, such as the deaminases or the phosphorylases. We set the flux of adenosine deaminase to 0.01 mM/hr as it a very low flux enzyme based on an earlier model (Joshi, 1990).
This assignment sets the flux of $\textbf{p}_{2}$ to zero and $\textbf{p}_{3}$ to 0.01 mM/hr. We pick the flux through $\textbf{p}_{2}$ to be zero since it overlaps with $\textbf{p}_{5}$ and gives flux values to all the reactions in the pathways.
With these pathays and numerical values, the steady state flux vector can be computed as the weighted sum of the corresponding basis vectors. The steady state flux vector is computed as an inner product:
# Set independent fluxes to determine steady state flux vector
independent_fluxes = {
ampsn.reactions.SK_amp_c: 0.0,
ampsn.reactions.SK_ade_c: -0.014,
ampsn.reactions.ADNK1: 0.12,
ampsn.reactions.SK_adn_c: -0.01,
ampsn.reactions.AMPDA: 0.014}
# Compute steady state fluxes
ssfluxes = ampsn.compute_steady_state_fluxes(
minspan_paths,
independent_fluxes,
update_reactions=True)
table_12_8 = pd.DataFrame(list(ssfluxes.values()), index=reaction_ids,
columns=[r"$\textbf{v}_{\mathrm{stst}}$"])
Table 12.8: The steady state flux through the AMP metabolism pathway.
and can be visualized as a bar chart:
fig_12_3, ax = plt.subplots(nrows=1, ncols=1, figsize=(10, 5))
# Define indicies for bar chart
indicies = np.arange(len(reaction_ids))+0.5
# Define colors to use
c = plt.cm.coolwarm(np.linspace(0, 1, len(reaction_ids)))
# Plot bar chart
ax.bar(indicies, list(ssfluxes.values()), width=0.8, color=c);
ax.set_xlim([0, len(reaction_ids)]);
# Set labels and adjust ticks
ax.set_xticks(indicies);
ax.set_xticklabels(reaction_ids, rotation="vertical");
ax.set_ylabel("Fluxes (mM/hr)", L_FONT);
ax.set_title("Steady State Fluxes", L_FONT);
# Add a dashed line at 0
ax.plot([0, len(reaction_ids)], [0, 0], "k--");
Figure 12.3: Bar chart of the steady-state fluxes.
We can perform a numerical check make sure that we have a steady state flux vector by performing the multiplication $\textbf{Sv}_{\mathrm{stst}}$ that should yield zero.
A numerical QC/QA: Ensure $\textbf{Sv}_{\mathrm{stst}} = 0$
pd.DataFrame(
ampsn.S.dot(np.array(list(ssfluxes.values()))),
index=metabolite_ids,
columns=[r"$\textbf{Sv}_{\mathrm{stst}}$"],
dtype=np.int64).T
The approximate steady state values of the metabolites are given above. The mass action ratios can be computed from these steady state concentrations. We can also compute the forward rate constant for the reactions:
percs = ampsn.calculate_PERCs(update_reactions=True)
Table 12.9: AMP Salvage Network enzymes, loads, transport rates, and their abbreviations. For irreversible reactions, the numerical value for the equilibrium constants is $\infty$, which, for practical reasons, can be set to a finite value.
These estimates for the numerical values for the PERCs are shown in Table 12.9. These numerical values, along with the elementary form of the rate laws, complete the definition of the dynamic mass balances that can now be simulated. The steady state is specified in Table 12.8.
The AMP degradation and biosynthetic pathways determine the amount of AMP and its degradation products present. We thus define the 'AMP charge,' calculated as the AMP concentration divided by the concentration of AMP and its degradation. This pool consists of IMP, ADO, INO, R1P, R5P, and PRPP. The AMP charge is
$$\begin{equation*} r_{\text{AMP}} = \frac{\text{AMP}}{\text{AMP} + \text{IMP} + \text{ADO} + \text{INO} + \text{R1P} + \text{R5P} + \text{PRPP}} \end{equation*}$$The AMP charge, $r_{AMP}$, is a little over 50% at steady-state calculated below. Thus, about half of the pentose in the system is in the AMP molecule and the other half is in the degradation or biosynthesis products.
The maintenance of AMP at a physiologically meaningful value is determined by the balance of its degradative and biosynthetic pathways. In Chapter 10 we saw that a disturbance in the ATP load on glycolysis leads to AMP exiting the glycolytic system, requiring degradation. Here we are interested in seeing how the balance of the AMP pathways is influenced by a sudden change in the level of AMP.
We thus simulate the response to a sudden 10% increase in AMP concentration as a proxy for a disturbance that increases AMP concentration, such as the simulations performed in the previous chapters.
t0, tf = (0, 1000)
sim_ampsn = Simulation(ampsn)
conc_sol, flux_sol = sim_ampsn.simulate(
ampsn, time=(t0, tf, tf*10 + 1),
perturbations={"amp_b": "amp_b * 1.1"})
conc_sol.make_aggregate_solution(
"AMP_Charge_Ratio",
equation="(amp_c) / (amp_c + adn_c + ade_c + imp_c + ins_c + prpp_c + r1p_c + r5p_c)")
fig_12_4, axes = plt.subplots(nrows=1, ncols=2, figsize=(12, 4),
constrained_layout=True)
(ax1, ax2) = axes.flatten()
plot_time_profile(
conc_sol, observable=ampsn.metabolites, ax=ax1,
plot_function="loglog",
xlabel="Time [hr]", ylabel="Concentration [mM]",
title=("(a)Concentrations", L_FONT));
plot_time_profile(
conc_sol, observable="AMP_Charge_Ratio", ax=ax2,
plot_function="semilogx",
xlabel="Time [hr]", ylabel="Ratio", ylim=(0.4, .6),
title=("(b) AMP Charge Ratio", L_FONT));
Figure 12.4: (a) The concentrations of the AMP metabolism pathway and (b) The AMP charge after a sudden 10% increase in the AMP concentration at $t=0$.
The response of the degradation pathways to AMP increase is a rapid conversion of AMP to IMP followed by a slow conversion of IMP to INO that is then rapidly exchanged with plasma, Figure 12.5. This degradation route is kinetically preferred. A sudden decrease in AMP has similar, but opposite responses. If AMP is suddenly reduced in concentration, the flux through this pathway drops.
Figure 12.5: The dynamic response of the AMP synthesis/degradation sub-network to a sudden 10% increase in the AMP concentration at $t=0$. The most notable set of fluxes are shown.
The AMP metabolic sub-network of Figure 12.1 can be integrated with the glycolysis and pentose pathway network from the last chapter. The result is shown in Figure 12.6. This integrated network represents the core of the metabolic network in the red blood cell, the simplest cell in the human body.
Figure 12.6: The core metabolic network in the human red blood cell comprised of glycolysis, the pentose pathway, and adenine nucleotide metabolism. Some of the integration issues discussed in the text are highlighted with dashed ovals.
Given the many points of contact created between the AMP sub-network and the combined glycolytic and pentose pathway network, there are a few interesting integration issues. They are highlighted with dashed ovals in Figure 12.6 and are as follows:
In the sub-network described above, the stoichiometry of the PRPP synthase reaction is $$\begin{equation} \text{R5P} + 2\text{ATP} \rightarrow \text{PRPP} + 2\text{ADP} + \text{H} \tag{12.2} \end{equation}$$
but in actuality it is $$\begin{equation} \text{R5P} + \text{ATP} \rightarrow \text{PRPP} + \text{AMP} + \text{H} \tag{12.3} \end{equation}$$
this difference disappears since the ApK reaction is in the combined glycolytic and pentose pathway network:
i.e., if Equations (12.3) and (12.4) are added, one gets Eq. (12.2).
The ATP cost of driving the biosynthetic pathways to AMP is now a part of the ATP load in the integrated model.
The AMP exchange reaction disappears, and instead we now have exchange reactions for ADO, ADE, INO, and HYP, and the deamination reactions create an exchange flux for $\text{NH}_3$.
Just as in the last chapter, we start by loading all the models, and then combine them:
glycolysis = load_json_model(path.realpath(path.join(
"models", "SB2_Glycolysis.json")))
ppp = load_json_model(path.realpath(path.join(
"models", "SB2_PentosePhosphatePathway.json")))
ampsn = load_json_model(path.realpath(path.join(
"models", "SB2_AMPSalvageNetwork.json")))
fullppp = glycolysis.merge(ppp, inplace=False)
fullppp.id = "Full_PPP"
fullppp.remove_reactions([
r for r in fullppp.boundary
if r.id in ["SK_g6p_c", "DM_f6p_c", "DM_g3p_c", "DM_r5p_c"]])
fullppp.remove_boundary_conditions(["g6p_b", "f6p_b", "g3p_b", "r5p_b"])
core_network = fullppp.merge(ampsn, inplace=False)
core_network.id = "Core_Model"
Then a few obsolete exchange reactions have to be removed.
for boundary in core_network.boundary:
print(boundary)
core_network.remove_reactions([
r for r in core_network.boundary
if r.id in ["DM_amp_c", "SK_amp_c"]])
core_network.remove_boundary_conditions(["amp_b"])
We then need to change the PRPP synthesis reaction to have the correct stoichiometry and PERC:
# Note that reactants have negative coefficients and products have positive coefficients
core_network.reactions.PRPPS.subtract_metabolites({
core_network.metabolites.atp_c: -1,
core_network.metabolites.adp_c: 2})
core_network.reactions.PRPPS.add_metabolites({
core_network.metabolites.amp_c: 1})
core_network.reactions.PRPPS.kf = 0.619106
The merged model contains 40 metabolites and 45 reactions:
print(core_network.S.shape)
Note that the merged model is not in a steady-state:
t0, tf = (0, 1e3)
sim_core = Simulation(core_network)
conc_sol, flux_sol = sim_core.simulate(
core_network, time=(t0, tf, tf*10 + 1))
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(12, 4))
(ax1, ax2) = axes.flatten()
plot_time_profile(
conc_sol, ax=ax1, legend="lower outside",
plot_function="loglog", xlabel="Time [hr]",
ylabel="Concentration [mM]", title=("Concentrations", L_FONT));
plot_time_profile(
flux_sol, ax=ax2, legend="lower outside",
plot_function="semilogx", xlabel="Time [hr]",
ylabel="Flux [mM/hr]", title=("Fluxes", L_FONT));
As in the previous chapter, we can perform a set of transformations to group the species and reactions into organized groups.
# Define new order for reactions
new_metabolite_order = [
"glc__D_c", "g6p_c", "f6p_c", "fdp_c", "dhap_c",
"g3p_c", "_13dpg_c", "_3pg_c", "_2pg_c", "pep_c",
"pyr_c", "lac__L_c", "_6pgl_c", "_6pgc_c", "ru5p__D_c",
"xu5p__D_c", "r5p_c", "s7p_c", "e4p_c", "ade_c", "adn_c",
"imp_c", "ins_c", "hxan_c", "r1p_c", "prpp_c", "nad_c",
"nadh_c", "amp_c", "adp_c", "atp_c", "nadp_c", "nadph_c",
"gthrd_c", "gthox_c", "pi_c", "h_c", "h2o_c", "co2_c", "nh3_c"]
if len(core_network.metabolites) == len(new_metabolite_order):
core_network.metabolites = DictList(core_network.metabolites.get_by_any(new_metabolite_order))
# Define new order for metabolites
new_reaction_order = [
"HEX1", "PGI", "PFK", "FBA", "TPI", "GAPD", "PGK", "PGM",
"ENO", "PYK", "LDH_L", "G6PDH2r", "PGL", "GND", "RPE",
"RPI", "TKT1", "TKT2", "TALA", "ADNK1", "NTD7", "ADA","AMPDA",
"NTD11", "PUNP5", "PPM", "PRPPS", "ADPT", "ADK1",
"ATPM", "DM_nadh", "GTHOr", "GSHR", "SK_glc__D_c", "SK_pyr_c",
"SK_lac__L_c", "SK_ade_c", "SK_adn_c", "SK_ins_c", "SK_hxan_c",
"SK_pi_c", "SK_h_c", "SK_h2o_c", "SK_co2_c", "SK_nh3_c"]
if len(core_network.reactions) == len(new_reaction_order):
core_network.reactions = DictList(core_network.reactions.get_by_any(new_reaction_order))
The stoichiometric matrix for the merged core network model is shown in Table 11.11. This matrix has dimensions of 40x45 and its rank is 37. The null is of dimension 8 (=45-37) and the left null space is of dimension 3 (=40-37). The matrix is elementally balanced.
We have used colors in Table 12.10 to illustrate the organized structure of the matrix.
Table 12.10: The stoichiometric matrix for the merged core network model in Figure 12.6. The matrix is partitioned to show the glycolytic reactions (yellow) separate from the pentose phosphate pathway (light blue) and the AMP salvage network (light green). The cofactors (light orange) and inorganics (pink) are also grouped and shown. The connectivities, $\rho_i$ (red), for a compound, and the participation number, $\pi_j$ (cyan), for a reaction are shown. The second block in the table is the product $\textbf{ES}$ (blue) to evaluate elemental balancing status of the reactions. All exchange reactions have a participation number of unity and are thus not elementally balanced. The last block in the table has the 8 pathway vectors (purple) for the merged model. These vectors are graphically shown in Figure 12.7. Furthest to the right, we display the time invariant pools (green) that span the left null space.